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1. Introduction 

The most celebrated stochastic differential equation in Physics is undoubtedly the 
Langevin equation 



describing (overdamped) Brownian motion, where r)(t) is a Gaussian white noise. 
Langevin equations are usually written down according to the following phenomeno- 
logical scheme: one starts with a deterministic equation, e.g. dx/dt = in the absence 
of any external force, and adds a noise term which mimics the interactions of the system 
with its external environment. The resulting Langevin equation is in general much more 
involved than the original deterministic equation. For a wide class of problems, Langevin 
equations however remain analytically tractable [1, 2, 3, 4, 5, 6]. 

Another broad class of dynamical systems of interest comprises memory. Pheno- 
mena involving memory effects and time delays are indeed ubiquitous. To be more 
specific, differential equations with delay play an important role in many fields of 
the physical sciences, technology, economy, and so on, whenever memory effects 



dx(t) 
dt 



T](t) 
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must be taken into account. They have been the subject of extensive mathematical 
studies [7, 8, 9]. In such a circumstance, one must know the whole past of the system, 
in addition to its present, in order to predict its immediate future. The simplest example 
of a differential equations with delay is provided by the linear equation 

^ = *«-!), (1.2) 

with a constant delay equal to 1. To solve this equation, we need to know the initial 
values x(t) on the interval [0,1], say. We can then recursively derive x(t) at all subsequent 
times. 

The goal of the present work is to investigate the joint effects of stochasticity and 
of memory by considering differential equations with unbounded random time delay. 
Because of this memory effect, the resulting stochastic dynamics are non-Markovian. 
The prototype of such an equation is the following first-order linear equation: 

^ = s(r(f)), (1-3) 

where r(t) is an earlier time (i.e., < r(t) < t), chosen according to some random 
process, the delay itself being the time difference t — r(t). At variance with the Langevin 
equation (1.1), no explicit noise appears in (1.3), so that the random delay process 
{r(t)} is the only source of stochasticity. In the following we focus our attention onto 
the simplest situation of a uniform sampling where, at each time t, the value of r(t) is 
chosen uniformly over the whole past, i.e., in the interval [0,t], independently of what 
occurs at other times. 

Discrete avatars of the present problem, namely linear recursions with unbounded 
discrete random delays, have a long history dating back to the pioneering investigation 
of the random Fibonacci sequence by Mark Kac [10]. Stochastic differential equations 
with unbounded random delay of the form (1.3) have however not been investigated 
so far, to the best of our knowledge. Several mathematical works have established a 
range of results on difference or differential equations with random delay [11, 12, 13, 14], 
concerning the stability and the ergodic behavior of such stochastic dynamical systems. 
More specifically, the case of differential equations with a stationary random delay is 
dealt with in [13], the emphasis being on the existence of a stationary random solution 
and on its stability under discretization, whereas a multiplicative ergodic theorem is 
proved in [14] for a class of difference equations with random delay. The relationship 
between the latter investigations and the present work is however tenuous. It is worth 
emphasizing that the random delay r(t) considered here is not a stationary process, 
because it is uniformly distributed in the growing time interval [0, t}. The general results 
derived in [11, 12, 13, 14] therefore do not apply in the present case. In particular, 
we shall find that the typical growth of the function x(t) solving (1.3) is a stretched 
exponential of the form exp(2y / t), whereas it is proved in [14] that for a stationary 
random delay the growth is always exponential. Finally, yet another field of application 
of random delays can be found in recent generalizations of the Black-Scholes model, 
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where the non- uniformity of the information in the market and the time delay until new 
information becomes generally available are taken into account [15]. 

The setup of this article is as follows. Section 2 contains a detailed study of the 
first-order linear equation (1.3). Our main result is that the solution to this equation is 
self-averaging. Equation (1.3) will indeed be shown to be equivalent to the deterministic 
integro-differential equation (2.18), which boils down to the second-order differential 
equation (2.19). In Section 3 we investigate a few other examples of dynamical systems 
with random delay, involving either higher-order time derivatives or non-linearities. 
Section 4 is devoted to the discrete counterpart of (1.3), namely the random linear 
recursion 

defining a random Fibonacci sequence. The special case where e — 1, originally 
considered by Kac [10], has also been studied in [16], whereas the case of an arbitrary 
time step e has been investigated in [17]. Here we put the main emphasis on the crossover 
between the fluctuating discrete problem and the deterministic continuous one as the 
time step e goes to zero. Section 5 contains a brief discussion of our findings. 



2. A first-order equation with random delay 

In this section we investigate the first-order differential equation (1.3) in full detail. As 
the latter equation is linear, it suffices to consider the initial condition x(0) = 1. 



2.1. Examples of deterministic delay 

In order to forge our intuition, let us first consider a few examples where the delay 
process r(t) is deterministic. 

• In the absence of delay, i.e., r(t) = t, we have the ordinary differential equation 
dx/dt = x, whose solution is a pure exponential: 

x(t) = e*. (2.1) 

• With a constant delay a, i.e., r(t) = t — a, we obtain an equation similar to (1.2). 
The initial data must be the full function x(t) for < t < a. The solution still grows 
asymptotically exponentially, albeit with a reduced rate, as 

x{t) ~ e w *, (2.2) 

where oj < 1 is the real solution of the characteristic equation 

u = e-^. (2.3) 

The rate u decreases monotonically as a function of the delay a. For a C 1, we have 
uj = 1 — a + 3a 2 /2 + • • • In the opposite regime (a ^> 1), we have u> ~ (In a) /a. The 
characteristic equation (2.3) also has an infinite sequence of complex roots, describing 
the damped oscillations displayed by x(t) for generic initial data, before the asymptotic 
exponential growth sets in. 
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• For an arbitrarily slowly varying delay a(t) = t — r(t), the solution can be argued to 
grow as 

x(t) ~ exp ^jf* u{s) ds^j , (2.4) 

where the instantaneous growth rate u(t) is related to the delay a(t) as 

u (t) = e - u ®« t \ (2.5) 

so that oj(t) ~ (\n a(t)) / a(t) is small whenever a(t) is large but slowly varying. 

• For a delay growing linearly in time, i.e., r(t) = bt, where b < 1 is fixed, the solution 
to (1.3) has the exact series representation 

x(f) = £&* ( *- 1)/2 ^. (2.6) 
k>o k. 

The saddle-point method yields the asymptotic growth law 

*® ~ exp (S) ' (27) 

which is faster than any power law, but slower than any stretched exponential. 

The growth law in all these cases is vastly different, so it is not clear what to expect 
in the stochastic case of a uniform random delay. We shall investigate the properties 
of equation (1.3) by successively calculating the average of the variable x(t) and its 
fluctuations. 

2.2. First moment 

Let us begin our analysis of (1.3) with a uniform random delay by studying the first 
moment (average) (x(t)). To do so, we are led to introduce two quantities: 

M(t) = (x(t)), N(t)= [\x(T))dr. (2.8) 

Jo 

Here and throughout the following, brackets (• • •) denote an average over the realizations 
of the random delay process {r(t)}. The above quantities obey the differential equations 

dM N dN , nns 

— = — , —r— — M, 2.9 
dt t dt v ' 

with initial values M(0) = 1, JV(0) = 0. 

The average M(t) therefore obeys the second-order differential equation 

d / dM\ d 2 M dM 

— t—— = t — —j— H — — = M. 2.10 
dt\ dt J dt 2 dt y ' 

The latter equation can be solved explicitly by looking for the series expansion of M(t). 
Its solution is 

-i-k 

M(t) = J2 — = I (2Vi) : (2.11) 

fc>0 K - 

where I is the modified Bessel function. We note that equation (2.10) can alternatively 
be transformed into a Bessel equation by using the change of variable T = 2y/i. 
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The average solution exhibits a stretched exponential growth of the form 

M(t) « j= - - . (2.12) 



Putting this growth law in perspective with (2.4) and (2.5), we are led to conclude that 
the relevant values of the time delay lie in the range a(t) — t — r(t) ~ \ft. 

2.3. Second moment 

In order to explore the distribution of x(t), we now turn to the second moment (x(t) 2 ). 
In analogy with (2.8), we are led to introduce three quantities: 



Q(t) = (x(t) 2 }, R(t)= f (x(r)x(t)) dr, 

Jo 



t r t 



which obey 



S(t)= (x(r)x(r'))drdr', (2.13) 
Jo Jo 

with initial values Q(0) = 1, R(0) = S(0) = 0. 

It can easily be checked that the solution to the above equations is 

Q(t) = M(t) 2 , R(t) = M(t)N(t), S(t) = N(t) 2 . (2.15) 

Indeed, as a consequence of (2.9) and (2.14), both sides of each identity of (2.15) obey 
the same first-order differential equations, with the same initial values. The meaning of 
the above identities will be discussed in the next section. 



2.4- Deterministic behavior 

The first of the identities (2.15) tells us that we have identically (x(t) 2 } = (x(t)} 2 for all 
times t. This implies that the solution to the differential equation (1.3) with random 
delay is with certainty equal to the deterministic function 

x(t) = M(t) = I (2Vt). (2.16) 

This solution is compared in Figure 1 with the exponential solution (2.1) in the absence 
of delay. 

The solution (2.16) to the differential equation (1.3) is therefore self-averaging. This 
absence of fluctuations is striking at first sight. The stochastic differential equation (1.3) 
is indeed driven by the random delay process {r(£)}. 

A first observation to be made is that this self-averaging property is a peculiarity 
of the continuous-time limit. The discrete counterpart of the differential equation (1.3), 
i.e., the random linear recursion (1.4), indeed exhibits the generic kind of fluctuations 
to be expected in a stochastic dynamical system, including non-trivial Lyapunov 
exponents [10, 16, 17]. The crossover between the fluctuating discrete problem and 
the deterministic continuous one as the time step e goes to zero will be investigated in 
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Figure 1. Logarithmic plot of the solution x(t) to the differential equation (f.3). 
Black: no delay (see (2.1)). Red: random delay (see (2.16)). 



detail in Section 4. More generally, applying any discretization scheme to the stochastic 
differential equation (1.3) (to be necessarily used e.g. in a numerical analysis) will break 
the self- averaging property and resurrect statistical fluctuations. 

The peculiar self-averaging property of the continuous problem can be understood 
in two complementary ways. 

On the physical side, the stochastic differential equation (1.3) samples the process 
{r(t)} for infinitely many different times, however small t is. Fluctuations are 
therefore cut down by an infinite noise reduction factor, and so the resulting process 
is deterministic. In the presence of a discrete time step e, the random delay process 
at time t is only sampled a large but finite number of times, of order t/e, so that the 
reduced variance of the process x(t) is expected to be proportional to e. This prediction 
will be corroborated by the analysis of Section 4.2. 

On the mathematical side, our choice for the delay process, namely that r(i) is 
drawn at each time t independently of what occurs at all the other instants, implies 
that a typical realization of r(t) is very wild function of t, which is nowhere continuous, 
and not integrable. The formal short-time expansion of the solution, 



thus involves, as its first non-trivial term, an integral which is ill-defined. The most 
natural definition of the latter integral [18, 19] just consists in replacing it by its average 
value, i.e., t 2 /4, since (t(s)) = s/2. 

An efficient way of expressing the self-averaging property of the continuous problem 
is to say that (1.3) is equivalent to the deterministic integro-differential equation 




(2.17) 
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which boils down to the second-order differential equation 

d 2 x dx , n , _ 

( dF + dF^' (2 ' 19) 

To sum up, the self-averaging property is intrinsic to the continuous-time random 
delay process. It is therefore expected to hold quite generally for dynamical systems 
described by differential equations with random delay, involving either higher-order time 
derivatives or non-linearities. The net effect of the random delay is to increase the order 
of the differential equation by one unit. Several examples will be dealt with in Section 3. 

3. A panorama of equations with random delay 

3.1. Oscillations induced by delay in a first-order linear equation 

We start our panorama of differential equations with random delay by considering the 
equation 

^ = -x(r(t)), (3.1) 

which is obtained from (1.3) by changing the sign of the rate of evolution. This formally 
amounts to changing the sign of time. 

In the absence of delay, we have the differential equation dx/dt = —x, whose 
solution relaxes exponentially fast to zero, as 

x(t) = e~*. (3.2) 

With random delay, x(t) obeys the integro-differential equation 

dx(t) 1 rt 
dt t Jo 

It is therefore ruled by the second-order differential equation 

d 2 x dx 

«dP + a = -*- (34) 

whose solution is 

-t) k r cos(2 v / t - tt/4) 



x(t)6t. (3.3) 



*W = E^ = ^)^=^, (3.5) 



where Jo is the Bessel function. 

The solution x(t) exhibits a very slow power-law decay in 1/t 1 ^ 4 , modulated by 
oscillations which slow down in the course of time. In particular, x(t) vanishes at an 
infinite sequence of instants t k = j'f/4 w (n 2 /A)k 2 , where the j k are the zeros of J . We 
have h = 1.445796 . . ., t 2 = 7.617815 . . ., t 3 = 18.721751 . . ., and so on. 

Figure 2 shows a comparison between the solution (3.2) in the absence of delay, 
falling off exponentially fast to zero, and (3.5) with random delay, relaxing in a very 
slow and non-monotonic way. 
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Figure 2. Solution x(i) to the differential equation (3.1). Black: no delay (sec (3.2)). 
Red: random delay (see (3.5)). Blue: stable fixed point (x = 0). 



3.2. The harmonic oscillator with random delay 

The harmonic oscillator with random delay provides another illustration of the above 
framework. It is defined by the dynamical equation 

^ = -x{r{t)), (3.6) 

in reduced units. We assume the initial conditions are x(0) = 1, dx(0)/dt = 0. 
In the absence of delay, we have 

x(t) = cost. (3.7) 

The trajectory remains bounded. Its total energy, 

*(«) = 5 (»(«>' + < 3 - 8 > 

is conserved and equals E = 1/2. 
With random delay, x(t) obeys 

d 2 x(t) 1 r* . . , 

- / x(t) dr, (3.9) 



dt 2 t Jo 

so that 

d 3 x d 2 x ,„ , v 

The solution to this third-order differential equation is 

/ (-2t 2 ) k k\ (I 1 t 2 \ , 

*<*) = g i p^i-=^(5. f -g). ("I) 

where F 2 is a hypergeometric series. 
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The asymptotic behavior of x(t) at late times is drastically affected by the random 
delay. Setting z = —t 2 in the series expansion of (3.11), we obtain the following estimate 
by means of the saddle-point method, including its absolute prefactor: 

iW ^273 eXP (I" /3 )' (3 ' 12) 
For a positive time t, we have z 1 ^ 3 = (— t 2 ) 1 / 3 = e ±27rl /3^ 2 /3_ Summing the contributions 
of both complex conjugate saddle points, we end up with the following stretched 
exponential growth for the position, modulated by stretched oscillations: 




(3.13) 



Let us keep defining the total energy of the delayed harmonic oscillator by (3.8). 
This quantity is asymptotically dominated by its first (i.e., potential) term. It therefore 
grows as E(t) « x(t) 2 /2. The kinetic energy is relatively suppressed by a factor of 
order t~ 2 / 3 . Finally, the relevant values of the time delay are in the range a(t) = 
t - r(t) ~ t 1 ' 3 . 

Figure 3 shows a comparison between the periodic solution (3.7) describing the 
harmonic oscillator in the absence of delay and the growing solution (3.11) with 
random delay. 




t 



Figure 3. Solution x(t) describing the harmonic oscillator. Black: no delay (see (3.7)). 
Red: random delay (see (3.11)). Blue: fixed point (x = 0). 



3.3. Higher-order linear equations 

We now turn to the higher-order analogues of (1.3), namely 

^=*M*)), (3-14) 

where m > 1 is an arbitrary integer. For definiteness we assume that the initial condition 
is x(0) = 1, whereas the first m — 1 derivatives vanish at t — 0. 
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In the absence of delay, we have the differential equation d m x/dt m = x, whose 
solution is 



-imk i m—1 t 

±— = - Y exp(e 2 ^ m t) » -. 
fc>0 {mk)\ m ^ t^oo m 



*(t) = E 7-TTT = i- E exp(e 2 ^t) « (3.15) 



This solution grows asymptotically exponentially in i, with unit rate, irrespectively of 
the order m. 

With random delay, x(t) obeys 

d m x(t) 1 ft N , 

' x(r)dr, (3.16) 



dt m t Jo 
so that 

d m+1 x d m x , „. 

t - — + - — =x. 3.17 

The solution to the latter equation is 



X 



(t)= E 



fc>0 



_|_ J_^j m k (mk)\ 



/ 1 1 2 m - 1 t m \ 

— oFm ; — — -r J, (3.18) 

Vm mm m m m+1 J 

where oFm is a generalized hypergeometric series. The saddle-point method yields a 

stretched exponential growth of the form 

x(t) ~ f(™-2)/(2(m+i)) exp (=£1 . (3. 19 ) 

Both the stretching index m/(m+l) and the exponent (m — 2)/(2(m + l)) of the power- 
law prefactor increase with the order m. The relevant values of the time delay are in 
the range a(t) = t - r(t) ~ tV(«H-i). 

5.^. .A non-linear first- order equation with quadratic coupling 

In order to illustrate the effect of random delay on non-linear dynamical systems, we 
focus our attention onto the case study of a first-order equation with a quadratic non- 
linearity with a positive or a negative coupling. 

Positive coupling. In the case of a positive coupling, we consider the differential 
equation 

= *(f (*))', (3.20) 

in reduced units. For definiteness we set again x(0) = 1. 

In the absence of delay, we obtain the equation dx/dt = x 2 , whose solution 

x(t) = (3.21) 

blows up in a finite time. 
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In the presence of a random delay, the self-averaging property implies that x(t) is 
again a deterministic function and that it obeys 
dx(t) _ 1 f\ 

so that 



« -tJo X(T)dT - (3 ' 22) 

+ (3.23) 

At variance with the previous examples of linear equations, the non-linear differential 
equation (3.23) cannot be solved by analytical means. In qualitative analogy with (3.21), 
the solution to (3.23) can however be expected to diverge in a finite time. A local analysis 
indeed shows that x(t) exhibits a quadratic divergence of the form: 

x(t) = -^i- + -il--^ + 14 < t - 2 i °) + ... (3.24) 
; t-tfo (t-t ) 2 5(t-t ) 25t 125t§ 

The only unknown in the above expansion is the divergence time to, which depends on 
the initial condition. For x(0) = 1 we obtain t ~ 3.140857. The effect of a random 
delay is therefore that the solution blows up later, but faster (see Figure 4, left). 

Negative coupling. The case of negative coupling, i.e., 

^ = -x(r(t))\ (3.25) 

is formally obtained from the previous one by changing the sign of time. For definiteness 
we set again x(0) = 1. 

In the absence of delay, the solution 

x(t) = (3.26) 

falls off slowly to zero. 

In the presence of a random delay, x(t) obeys 

tfx dx = _^ 

dt 2 dt v 1 

The solution x(t) decreases monotonically from x(0) = 1, crosses zero at a finite time 
ti ~ 2.133528, and keeps decreasing until it diverges to — oo at a later finite time 
t ~ 17.00447. The behavior of x(t) as t — > t is still given by the expansion (3.24), up 
to a simultaneous sign change in t and to- The effect of random delay is more drastic in 
this case. Instead of relaxing to zero, the solution overshoots and diverges to — oo in a 
finite time. 

Figure 4 shows a comparison between the solutions in the absence of delay and with 
random delay, for a positive (left) and a negative (right) quadratic coupling. 




Figure 4. Solution of the first-order differential equation with quadratic coupling. 
Left: positive coupling (see (3.20)). Right: negative coupling (see (3.25)). Black: no 
delay (see (3.21), (3.26)). Red: random delay (see (3.23), (3.27)). Blue: asymptotes. 



4. The discrete problem of a random linear recursion 

This section is devoted to a detailed study of the crossover between the fluctuating 
discrete problem and the deterministic continuous one as the scale of discretization 
goes to zero. Introducing a discrete time step e in (1.3), we obtain the random linear 
recursion (see (1.4)) 

■Z-n+l X n -\- £X m , (4.1) 

leading to a random Fibonacci sequence. At every discrete time n (such that t = ne), 
the label m is drawn at random, uniformly among the n + 1 integers 0, . . . , n. Setting 
xq — 1 for definiteness, we have x\ = 1 + £, whereas x<i takes the values 1 + 2e and 
1 + 2e + e 2 with equal probabilities, and so on. 

As recalled in the Introduction, the above problem of random Fibonacci sequences 
has a long history, at least in the case e — 1. The latter was indeed investigated by Mark 
Kac, following a discussion with Stan Ulam. Kac obtained the result j3 2 = ^2(5 + vT7) 
for the growth rate of the second moment (see (4.24)), and later referred to it as a 
tremendous formula with a square root of 17 in it [10]. The case e — 1 has been further 
studied in [16], whereas the problem with an arbitrary time step e has been investigated 
in [17]. Our choice of (4.1) as a discretization of (1.3) follows all these earlier works. 

As already announced in Section 2.4, there is a qualitative difference between 
the continuous-time differential equation (1.3) and the discrete equation (4.1). The 
continuous problem has a self-averaging and deterministic solution, whereas the discrete 
one exhibits the generic kind of fluctuations to be expected in a stochastic dynamical 
system. In the following we investigate many facets of the crossover between the 
fluctuating discrete problem and the deterministic continuous one. 

The following three regimes are worth being considered: 
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• Regime I: e — >■ at fixed n (and so t — > 0), 

• Regime II: n — )■ oo at fixed £ (and so t — > oo), 



• Regime III: £ — > and n — > oo at fixed t = ne. 

Let us anticipate the following stretched exponential growth of the moments in 
Regime II: 

(x k n ) ~ & = 2k lk ^ (4.2) 

The reduced growth rates 7^ are referred to as the (generalized) Lyapunov exponents. 

4-1. First moment 

In analogy with the analysis of the continuous problem performed in Section 2, we start 
by studying the first moment (x n ). To do so, we introduce the quantities 

n 

M n = (x n ), N n =J2(xrn), (4.3) 



m=0 



which obey the recursions 



M n+1 = M n + — ^—N n , N n+1 = M n+1 + N n , (4.4) 
n + 1 

with M = N = 1, hence 

(n + l)M n+1 - (2n + 1 + e)M n + nM n ^ = 0. (4.5) 

Equations (4.3)-(4.5) are the discrete analogues of (2.8)-(2.10). They can be solved 
recursively. We thus obtain 

e 2 3e 2 e 3 

Mi = l + e, M 2 = l + 2e+-, M 3 = l + 3e + — + —, (4.6) 

and so on, so that M n is a polynomial of degree n in e. 

Quantitative results can be derived in the three regimes defined above. 

• In Regime I, the M n have a regular expansion in the time step e: 

M n = l + ne+ n ±^ + ... (4.7) 

• In Regime II, it is convenient to introduce the generating series 

F(z) = M n z n . (4.8) 

n>0 

Equation (4.5) is equivalent to the differential equation 

dF 

(l-z) 2 ^ = (l-z + e)F, (4.9) 
whose properly normalized solution is 

p ez/{l-z) 

F(z) = ^— — . (4.10) 

This expression can be recognized as being the generating series of the Laguerre 
polynomials [20]. We thus recover the result M n = L n (—e) [17]. 
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The asymptotic behavior of the contour integral representation 

Mn =l—^- / (4.11) 

can be estimated by means of the saddle-point method. We thus obtain 

„-e/2 „2^/n£ 

M " K 7C(^- (4 ' 12) 



This expression coincides with the asymptotic behavior of the solution M(t) of the 
continuous problem (see (2.12)), up to an absolute multiplicative prefactor e~ £//2 . With 
the notation of (4.2), the Lyapunov exponent of order one therefore reads identically 

7i = 1- (4.13) 

The leading stretched exponential behavior of (4.12) can be recovered without 
solving the difference equation (4.5) exactly. A more direct route can indeed be taken 
if a stretched exponential behavior of the form M n ps A n e ftv/ " is assumed, where A n 
stands for a slowly varying prefactor. Equations (4.4) imply that a consistent hypothesis 
is N n ps Bny/ne^ 1 ^. Taking the continuum limit of (4.4), we obtain 

(jA n neB n , ^B n ziA n . (4.14) 
As a consequence, pi/2 is the eigenvalue of the 2x2 matrix 

*-(;;) (4.i5) 

whose real part is the largest. We thus recover 0i = 2-y/e. This line of thought will be 
used for higher moments, where an exact solution in terms of generating series will not 
be available. 

• In Regime III, the continuum limit of the difference equation (4.5) yields the 
differential equation (2.10), so that the M n converge to the solution M(t) (see (2.11)). 
This solution matches the behavior (4.7) in Regime I, as M(t) = 1 + t + t 2 /A + • • • as 
t — > 0. The behavior of M(t) as t — > oo also nicely matches (4.12), as already noticed. 

4-2. Second moment 

In order to study the second moment (x^), we have to introduce the four quantities: 

n 

S n = E (x,x m >, T n = Y,{x 2 m ). (4.16) 

The variance of x n therefore reads 

V n = v&rx n = Q n - Ml (4.17) 
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The above quantities obey the recursions 

2e e 2 e 

Qn+l = Qn H ; T-Rn H ; rT n , Rn+l = Qn+l + R-n H ; T^m 

n + l n + 1 n+1 

S n +1 = —Qn+l + 2-Rn+l + S n , T n+ \ = Q n+ i + T n , (4-18) 

with Q = R = S = T = I. 

Four quantities are needed to evaluate the second moment in the discrete problem, 
while three were sufficient in the continuous one (see (2.13)). The case of higher moments 
will be dealt with in Section 4.3. 

The above equations can be solved recursively. We thus obtain 

Qi = 1 + 2e + e 2 , Q 2 = 1 + As + 5e 2 + 2e 3 + —, 

Q 3 = 1 + 6e + 12e 2 + '-f- + ±f- + e 5 + £ -, (4.19) 

so that 

V 1 = 0, V 2 = -, V 3 = — + - + — , (4.20) 
' 4 ' s 12 2 36 V ; 

and so on. In general Q n and V n are polynomials of degree 2n in e. 

Many quantitative results, most of which are novel, can be derived in the three 

regimes defined above. 

• In Regime I, it can be shown that V n behaves as e 4 . Skipping details, let us give the 
result 

Vn ^-w^y\ ... (4 . 21) 

• In Regime II, the growth of the second moment can be studied by means of the 
approach sketched at the end of Section 4.1. Assuming a stretched exponential behavior 
of the form Q n ps C n e^ 2 ^, we are left after some algebra with the condition that l3 2 /2 
is an eigenvalue of the 4x4 matrix 

/0 2s e 2 \ 
1 e 
2 
\1 J 

The combination z = /3f/(4e) = A'j 2 obeys the quadratic equation 

z 2 - (A + e)z + 2e = 0, (4.23) 
so that the Lyapunov exponent of order two is 



(4.22) 



72 = ^ = iV 2 ( 4+£+v/T ^ T ^) • (4 - 24) 



This result was first obtained in the case e = 1 by Kac [10], and then for an arbitrary e 
in [17]. For small e, the Lyapunov exponent admits the expansion 

£ 3e 2 3e 3 , inrS 

72 = 1 + — + + • • • 4.25 

' 16 512 8192 y J 
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The reduced second moment therefore scales as 

( x n) ^ e 4(72-7i)v / ^e ^ (4.26) 



(Xn) 2 

for e <C 1 in Regime II (the e — >■ limit is taken after the n — > oo limit). 

• In Regime III, by taking the continuum limit of the recursions (4.18), we predict the 
following scaling behavior for the variance of the process: 

V n meW(t). (4.27) 

The variance therefore grows linearly with e, for all values of t — ne. This result is 
in agreement with the argument on the noise reduction factor exposed in Section 2.4. 
The reduced variance scales as 

varx n V n , w . , w . 

^ 5? = _A Md r (t ), X (i ) = ^, (4.28) 

where M(t) is given in (2.11). 

The scaling function W(t) is found to obey the third-order differential equation 

d 3 W „ d 2 W „ „ ,dW rt „, 2 d 2 n d0 ,„ , AnnS 



where 



0(t) = I J* M(s) 2 ds - Q J* M(s) ds 



= V mi (4 30) 

The differential equation (4.29) can be interpreted as follows: the function 0(t), 
measuring the variance of the temporal dispersion of the mean solution M(s) up to 
time t, acts as a source for the noise variance W(t) of the process at time t. 

The solution to (4.29) can be written as an explicit power series. Skipping algebraic 
details, we give the result 

W(t) = E^^-, (4.31) 

fc>3 

with 

= fc(fc 2 + 9fc-4) 3 _ 

4(/c + l)(2/c-l) 8 V k k h y ! 

where H n = J2m=i ^/m is the n-th harmonic number. 
For small times, we obtain the expansions 

„, M t 3 t 4 59t 5 v , , t 3 t 4 809t 5 , inn . 

Wit) = — + — + + • • • , It = + + • • • 4.33 

w 36 72 18 000 w 36 24 18 000 v ; 

The leading term matches the behavior V n ~ n 3 e 4 /36 in Regime I (see (4.21)). 

In the opposite regime of late times, the saddle-point method yields 

*M = T-TT + " (4 ' 34) 
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The reduced variance therefore grows as 



var x r 



:Vt 



'ne° 



W " 4 - 4 ' < 4 ' 35 > 
This power-law growth crosses over to the stretched exponential one (4.26) when the 
reduced variance becomes of order unity. This takes place at a late time, of order 
t ~ l/e 2 , i.e., n ~ l/e 3 . 

Figure 5 shows a plot of the scaling function X(t) describing the reduced variance 
of the process. 




Figure 5. The scaling function X(t) describing the reduced variance of the process. 



4-3. Higher moments 

Let us now turn to the higher moments (x^), where k is an arbitrary integer order. The 
approach exposed in the previous sections generalizes as follows. We are led to introduce 
a certain number of quantities which, together with (x^), obey closed linear recursion 
relations generalizing (4.4) and (4.18). The number of those quantities can be evaluated 
as follows. Each of them involves the product of x J n , for some j = 0, . . . , k, and of k — j 
factors x mi , where the labels rrii are not necessarily distinct. Attributing these labels 
amounts to partitioning the integer k — j. For fixed k and j, the number of distinct 
quantities is therefore Pk-j, where pk is the number of partitions of the integer k [20]. 

Hence the total number of quantities to be considered at order k is 

k 

ft = X>i- ( 4 - 36 ) 

Focusing our attention onto Regime II, and looking for a stretched exponential 
growth of the form (x^) ~ eP*v™ (see (4.2)), we are left with the condition that /3&/2 
is an eigenvalue of a square matrix of size x q k . The eigenvalues of £& occur in 
pairs of opposite numbers. The intuitive reason for that is that the coefficients of the 
recursion relations are rational in n, so that only even powers of y/n, i.e., even powers 
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k 


Pk 


qk 


A fc 





1 


1 


2 


1 




2 


2 


4 


2 




3 


3 


7 


3 


* 


4 


5 


12 


6 




5 


7 


19 


9 


* 


6 


11 


30 


15 




7 


15 


45 


22 


* 


8 


22 


67 


33 


* 


9 


30 


97 


48 





Table 1. Numbers involved in the analysis of the moment (x„) up to k = 9: pt is 
the number of partitions of the integer fc, q k is the number of quantities needed to 
obtain closed recursion relations, A k is the degree of the polynomial equation obeyed 
by the combination z = /3|/(4e) = (kjk) 2 - A star in the last column signals that the 
matrix H k has a zero eigenvalue. 



of /3k, should matter. If q^ is even, there are A& = q^/2 pairs of non-zero eigenvalues. 
If qt is odd, there are A& = (qk — l)/2 pairs of non-zero eigenvalues, whereas is a 
simple eigenvalue. As a consequence, the combination z = /3%/(4e) = (k'-fk) 2 obeys a 
polynomial equation of degree A&, generalizing (4.23). Table 1 gives the numbers pk, qk, 
and A k for the first few values of the order k. 
Let us now give some explicit results. 

• For k — 3, A 3 = 3, and the polynomial equation for z = [51/ (Ae) = 9^ is 

z 3 - (11 + Ae + e 2 )z 2 + (19 + 18e + lie 2 + e 3 )z 

- (9 - 18e + 16e 2 + 3e 3 ) = 0. (4.37) 

• For k = 4, A 4 = 6, and the polynomial equation for z = (31/ (Ae) = I67 2 is 

z 6 - (25 + 12e + he 2 + e 3 )z 5 

+ (168 + lQ8e + 12Ls 2 + 54e 3 + 10e 4 + e 5 )z 4 

- (400 + 492e + 487e 2 + 474e 3 + 174e 4 + 37e 5 + 3e 6 )z 3 

+ (256 + 1104e - 212e 2 + 7l2e 3 + 556e 4 + 251e 5 + AAs 6 + 2e 7 )z 2 

- 2e(384 - 22Ae - 92e 3 + 124e 4 + 55e 5 + Ae % )z 

+ 24e 4 (2 - e) 2 = 0. (4.38) 



For e = 1, besides f3 2 = ^2(5 + y/Vf) = 4.271558 . . ., we obtain /3 3 = 6.922223 . . . 
and (5 A = 10.120583... 

For small e, the Lyapunov exponents admit the expansion (see (4.25)) 

e 59e 2 1315 3 3e 9le 2 Al3e 3 

73 = 1+ 8 + 2304 + 36864 + '"' 74 = 1+ 16 + 1536 + 24576 + " 
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This behavior is fully analogous to what is observed in a broad class of disordered 
systems, prototypes of which are noisy dynamical systems [21] and the Anderson 
localization problem in one dimension [22, 23, 24]. The time step e plays the role 
of the strength of disorder, measured e.g. by the variance of the site energies in the case 
of the Anderson model with diagonal disorder. In analogy with the latter situation, 
we are tempted to deduce from the above results that the expansion of the Lyapunov 
exponent of arbitrary order k (not necessarily an integer) involves polynomials in k of 
increasing degrees, i.e., 

, (k - l)e (k - l)(32k - 37)e 2 

Ik = 1 + — + 1 (4.40) 

lk 16 4608 V ; 

In the opposite regime of large e, our explicit results for k — 2, 3, and 4 suggest 

the scaling behavior 

£ (fc-l)/2 

lk ~ — £ — • (4.41) 
4-4- Typical sequence and fundamental Lyapunov exponent 

The scaling form (4.2) of the moments implies that the mean logarithm of x n grows as 

(\nx n ) « 27^, (4.42) 

where 7 = 70 is the fundamental (usual) Lyapunov exponent. The latter describes the 
asymptotic growth law of the most probable or typical sequence: 

Mtyp ~ exp(27v/ne). (4.43) 

The exact calculation of the fundamental Lyapunov exponent 7 is beyond the reach of 
the present work. We can however write down its expansion at small e by setting k — 
in our conjectured formula (4.40). We thus obtain 

7=1 1 1 (4.44) 

1 16 4608 V ; 

Figure 6 shows a plot of numerical data for the Lyapunov exponent 7, obtained by 

means of a direct simulation of the random recursion (4.1). For each value of e, lnx n has 

been averaged over 10 7 different realizations of 10 4 steps each, and the outcome fitted 

to (4.42). For e = 1 our estimate 7 0.9448, hence 27 1.8896, fully corroborates the 

value 1.889 given in [16]. The data are in very good agreement with a fit incorporating 

the three terms of (4.44), as well as the falloff as e -1 / 2 suggested by (4.41). 



5. Discussion 

In this work we have investigated the joint effects of memory (i.e., time delay) and of 
extrinsic stochasticity (i.e., randomness) on the example of differential equations with 
unbounded random time delay. Our initial motivation was to investigate continuous- 
time analogues of the random Fibonacci sequences first considered by Kac [10]. 
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Figure 6. Numerical data for the Lyapunov exponent 7 against the time step e. Red 
line: rational fit incorporating the three terms of the expansion (4.44) and the falloff 
as e -1 / 2 suggested by (4.41). 

Stochastic differential equations with random time delay of this kind display a self- 
averaging property which leads to an unexpected deterministic behavior. A physical 
explanation of this striking self-averaging behavior, described in Section 2.4, relies on the 
fact that continuous time is infinitely divisible. So, during any time interval, whatever 
small, the random delay process {r(t)} will assume all possible values because of infinite 
sampling, and so the delay term will contribute only via its mean value. In other words, 
the averaging method used in classical mechanics [25] is exact for the present problem. A 
more mathematically inclined reader would argue that an equation with random delay of 
the type (1.3) is ill-defined because it identifies the derivative of a function (on the left- 
hand side) that always satisfies some smoothness properties (such as the intermediate 
value theorem) to a wild random function (on the right-hand side) with no smoothness 
property whatsoever; therefore, the random term needs to be regularized by replacing 
it by its almost-sure average value. 

In any case, randomly-delayed analogues of classical dynamical systems such as the 
first-order linear equation, the harmonic oscillator, or the non-linear growth model, have 
been shown to behave very differently from their deterministic counterparts. 

Our motivation to study continuous analogues of random Fibonacci sequences 
was to calculate the Lyapunov exponent that characterizes their typical growth by 
using differential methods (which are usually more versatile and powerful than discrete 
techniques). Fluctuations are however lost when taking the formal limit of a continuous 
time, so that a finite time step e has to be kept in order to preserve stochastic 
effects. We have investigated the various facets of the crossover between the fluctuating 
discrete problem and the deterministic continuous one. This led us, among many other 
outcomes, to the expansion (4.40) for the generalized Lyapunov exponents that measure 
the growth of the fc-th moment. The k — > limit of the latter result provides us 
with the expansion (4.44) for the fundamental Lyapunov exponent. Nevertheless, an 
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exact calculation of the latter quantity for random Fibonacci sequences still remains a 
challenging open problem. 

To close, it is worth mentioning that the potentially spectacular effects of long- 
ranged memory on random walks have been scrutinized in a long series of recent 
works [26, 27, 28, 29, 30, 31, 32]. 
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